Load network and scent data

From Tables S1 and S2 in: Kantsa, A., R.A. Raguso, A.G. Dyer, J.M. Olesen, T. Tscheulin, and T. Petanidou. 2018. Disentangling the role of floral sensory stimuli in pollination networks. Nature Communications 9: 1041. doi:[10.1038/s41467-018-03448-w](https://doi.org/10.1038/s41467-018-03448-w)

#load network
net <- read.delim("Kantsa_etal_Supp_Data1.csv", skip = 2)
rownames(net) <- net[,1]
net[,1] <- NULL
net <- as.matrix(t(net))
polls <- colnames(net)
snet <- sortweb(net) #sort by row & column totals

#load floral scents
scent.file <- "Kantsa_2018_Supp_Data2.xlsx"
plants <- getSheetNames(scent.file)
scentlist <- lapply(2:length(plants), read.xlsx, xlsxFile=scent.file)
plants <- plants[-1]
scent <- rbindlist(scentlist, fill=T, idcol="plant")[,-4]
colnames(scent)[2:3] <- c("cmpd","mean")
scent$plant <- factor(plants[scent$plant])
scent$cmpd <- factor(tolower(scent$cmpd))
cmpds <- levels(scent$cmpd)

scent.net <- dcast(scent, plant ~ cmpd, value.var="mean", fun.aggregate=mean, fill=0) #long to wide
scent.net <- as.data.frame(scent.net)
rownames(scent.net) <- scent.net[,1]
scent.net[,c(1,380)] <- NULL #weird NA column has one entry
scent.net <- as.matrix(scent.net)

scent.net.short <- scent.net[,colSums(scent.net)>0.1]

Plant phylogeny

#plants.class <- classification(plants, db="gbif")
#save(plants.class, file="plantsclass.Rdata")
load("plantsclass.Rdata")
plants.classtree <- class2tree(plants.class) #taxonomic tree
plants.tree <-phylomatic(taxa=names(plants.class), get = 'POST', storedtree="zanne2014")
plants.tree$tip.label[order(plants.tree$tip.label)] <- plants
par(mar=c(0,0,0,0))
cp <- cophylo(plants.tree, plants.classtree$phylo)
   Rotating nodes to optimize matching...
   Done.
plot(cp)

Pollinator phylogeny

polls <- gsub("Bug", "Hemiptera", polls, fixed=T)
polls <- gsub("Auchenorrhyncha", "Hemiptera", polls, fixed=T)
polls <- gsub("Fly", "Diptera", polls, fixed=T)
polls <- gsub("Beetle", "Coleoptera", polls, fixed=T)
#polls.class <- classification(polls, db="gbif")
#save(polls.class, file="pollsclass.Rdata")
load("pollsclass.Rdata")
polls.gbif <- factor(sapply(polls.class, function(x) tail(x$name, 1)))
net.lump <- aggregate(.~ polls.gbif, as.data.frame(t(net)), sum)
rownames(net.lump) <- net.lump[,1]
net.lump[,1] <- NULL
net.lump <- as.matrix(t(net.lump))
rownames(net.lump) <- sort(plants.tree$tip.label)
polls.classtree <- class2tree(unique(polls.class)) #taxonomic tree
polls.matched <- tnrs_match_names(lapply(polls.class, function(x) tail(x$name, 1)), context_name = "Animals")
   Warning: Some names were duplicated: 'andrena', 'asilidae', 'coleoptera',
   'hemiptera', 'buprestidae', 'buprestidae', 'calliphoridae', 'cerambycidae',
   'chrysididae', 'crabronidae', 'eumenidae', 'eumenidae', 'eumenidae',
   'diptera', 'diptera', 'diptera', 'diptera', 'diptera', 'hoplitis',
   'hoplitis', 'lasioglossum', 'noctuidae', 'oedemera', 'osmia', 'osmia'.
   Warning in check_tnrs(res): chrysura circe are not matched
polls.ids <- ott_id(polls.matched)
polls.ids <- polls.ids[!polls.ids %in% c(968359, 707236, 707875, 63881, 458953, 42699, 742352, 340559, 672497, 34898, 332969,356987)]
polls.tree <- tol_induced_subtree(ott_ids=polls.ids)
   
Progress [-----------------------------------] 0/7 (  0%) ?s
Progress [==================================] 7/7 (100%)  0s
                                                            
   Warning in collapse_singles(tr): Dropping singleton nodes with labels:
   Colletes ott115493, Crabronidae ott372234, Scolia ott881314, Coleoptera
   ott865243, Diptera ott661378, Paragus ott973598, Hemiptera ott603650
par(mar=c(0,0,0,0))
plot(polls.classtree, cex=0.5)

plot(polls.tree, cex=0.5)

Plant-pollinator network

Cophylogeny

mnet.lump <- melt(net.lump)
mnet.lump <- mnet.lump[mnet.lump$value>0,]
cophyloplot(plants.tree, polls.classtree$phylo, mnet.lump[,1:2], lwd=sqrt(mnet.lump[,3])/4, show.tip.label=F,  col='steelblue', length.line=0, gap=-20, space=60)

Pollinator groups

polls.cut <- dendextend::cutree(polls.classtree$phylo,5)
polls.order <- setNames(factor(sapply(unique(polls.class), function(x) x$name[4])), names(polls.class[!duplicated(polls.class)]))
polls.order.all <- setNames(factor(sapply(polls.class, function(x) x$name[4])),names(polls.class))
#table(polls.cut, polls.order)

net.grp <- aggregate(.~ polls.order.all, as.data.frame(t(net)), sum)
rownames(net.grp) <- net.grp[,1]
net.grp[,1] <- NULL
net.grp <- as.matrix(t(net.grp))
rownames(net.grp) <- sort(plants.tree$tip.label)

plants.majorpoll <- setNames(factor(colnames(net.grp)[apply(net.grp, 1, which.max)]), rownames(net.grp))

pcols <- setNames(brewer.pal(4,"Set1"), levels(plants.majorpoll))
plants.polltrees<-make.simmap(plants.tree,plants.majorpoll,model="ARD",nsim=100)
   make.simmap is sampling character histories conditioned on the transition matrix
   
   Q =
                 Coleoptera       Diptera  Hymenoptera Lepidoptera
   Coleoptera  -0.008667015  0.0006377689  0.008029247  0.00000000
   Diptera      0.017012461 -0.0224322446  0.005419784  0.00000000
   Hymenoptera  0.000000000  0.0000000000 -0.071181044  0.07118104
   Lepidoptera  0.078649350  0.0796271485  1.270049241 -1.42832574
   (estimated using likelihood);
   and (mean) root node prior probabilities
   pi =
    Coleoptera     Diptera Hymenoptera Lepidoptera 
          0.25        0.25        0.25        0.25
   Done.
plants.polltree<-summary(plants.polltrees,plot=FALSE)
plot(plants.polltree,colors=pcols,fsize=0.8,cex=c(0,0.3))
add.simmap.legend(colors=pcols,x=0.9*par()$usr[1],
    y=0.2*par()$usr[4],prompt=FALSE,fsize=0.9)

pcols.all <- setNames(c(pcols[1:2],brewer.pal(5,"Set1")[5],pcols[3:4]), levels(polls.order.all))

Heatmap

net.cut <- net[rowSums(net)>0,]
heatmaply(net.cut, scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
          scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
          distfun=vegdist, 
          hclustfun=function(x) hclust(x, method="mcquitty"), 
          row_side_colors = data.frame(P=plants.majorpoll[rowSums(net)>0]),
          row_side_palette=function(n) pcols, 
          col_side_colors=data.frame(PollinatorOrder=polls.order.all), 
          col_side_palette=function(n) pcols.all)
   Warning in heatmaply.heatmapr(hm, colors = colors, limits = limits,
   scale_fill_gradient_fun = scale_fill_gradient_fun, : The hover text for
   col_side_colors is currently not implemented (due to an issue in plotly).
   We hope this would get resolved in future releases.

Matrix plot

#interaction matrix plot
visweb(net, type="nested", circles=T, boxes=F, circle.max=1.2)

Web plot

#interaction web plot
plotweb(sqrt(net), empty=F, text.rot=90, method="cca", col.high=pcols.all[polls.order.all], col.low=pcols[plants.majorpoll], col.interaction=pcols.all[polls.order.all], bor.col.interaction=pcols.all[polls.order.all], bor.col.high=pcols.all[polls.order.all], bor.col.low=pcols[plants.majorpoll])

Network topography

networklevel(net)
                         connectance                     web asymmetry 
                          0.06312657                        0.63106796 
                   links per species            number of compartments 
                          1.95631068                        2.00000000 
               compartment diversity               cluster coefficient 
                          1.07901412                        0.02631579 
                          nestedness               weighted nestedness 
                          4.47259022                        0.58578362 
                       weighted NODF    interaction strength asymmetry 
                         10.80348739                        0.20372363 
            specialisation asymmetry                   linkage density 
                         -0.04424002                        3.57255189 
                weighted connectance                      Fisher alpha 
                          0.01734248                       89.61252715 
                   Shannon diversity              interaction evenness 
                          3.27629195                        0.37393976 
        Alatalo interaction evenness                                H2 
                          0.25285628                        0.57919894 
                number.of.species.HL              number.of.species.LL 
                        168.00000000                       38.00000000 
   mean.number.of.shared.partners.HL mean.number.of.shared.partners.LL 
                          0.34944397                        1.18776671 
              cluster.coefficient.HL            cluster.coefficient.LL 
                          0.34062967                        0.21733898 
     weighted.cluster.coefficient.HL   weighted.cluster.coefficient.LL 
                          0.72421977                        0.31277076 
                    niche.overlap.HL                  niche.overlap.LL 
                          0.12366249                        0.14889538 
                     togetherness.HL                   togetherness.LL 
                          0.03019814                        0.02302573 
                          C.score.HL                        C.score.LL 
                          0.75923950                        0.74100046 
                          V.ratio.HL                        V.ratio.LL 
                          3.23242952                       17.02718828 
                      discrepancy.HL                    discrepancy.LL 
                        248.00000000                      250.00000000 
                 extinction.slope.HL               extinction.slope.LL 
                          4.79325165                        1.62192500 
                       robustness.HL                     robustness.LL 
                          0.80798877                        0.61786266 
       functional.complementarity.HL     functional.complementarity.LL 
                       8044.31456516                     7651.23911634 
                partner.diversity.HL              partner.diversity.LL 
                          0.94103312                        1.28483306 
                       generality.HL                  vulnerability.LL 
                          2.67434342                        4.47076037

Scent network

Heatmap

heatmaply(sqrt(scent.net), scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
          scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
          distfun=vegdist, 
          hclustfun=function(x) hclust(x, method="mcquitty"), 
          row_side_colors = data.frame(P=plants.majorpoll),
          row_side_palette=function(n) pcols)

Web plot

#interaction web plot
plotweb(sqrt(scent.net), text.rot=90, method="cca", col.high=pal[2], bor.col.high=pal[2], col.low=pcols[plants.majorpoll], bor.col.low=pcols[plants.majorpoll], col.interaction=pal[5], bor.col.interaction=pal[5])

Plant-scent-pollinator network

NMDS

scent.mds <- metaMDS(sqrt(scent.net), autotransform = F, trace=F)
#scent.op <- ordiplot(scent.mds, type="n")
#points(scent.op, what="species", pch=3, col="grey50")
#points(scent.op, what="sites", cex=1.5, pch=19, col=pcols[plants.majorpoll])
#text(scent.op, what="sites", cex=0.8, col=pcols[plants.majorpoll])

Phylochemospace

scent.mds.taxa <- scent.mds$points
rownames(scent.mds.taxa) <- sort(plants.tree$tip.label)
rownames(scent.mds.taxa) <- gsub(" ", "_", rownames(scent.mds.taxa), fixed=T)
plants.tree$tip.label <- gsub(" ", "_", plants.tree$tip.label, fixed=T)

#plants.tree$phylo$edge.length <- plants.tree$phylo$edge.length+0.00001

#Do generalists cluster?
plants.numlinks  <- rowSums(decostand(net.lump, "pa"))
tip.pcols <- setNames(colorRampPalette(c("white", "darkblue"))(67)[plants.numlinks+1], 
                      sort(plants.tree$tip.label))
node.pcols<- setNames(c(tip.pcols[plants.tree$tip.label], rep("black",plants.tree$Nnode)),
                      1:(length(plants.tree$tip)+plants.tree$Nnode))

phylomorphospace(plants.tree, scent.mds.taxa, control=list(col.node=node.pcols), node.size=c(0,2), label="off", lwd=1)
text(scent.mds.taxa[,1], scent.mds.taxa[,2], rownames(scent.mds.taxa), cex=0.8, offset=0.5, xpd=T, col=tip.pcols)

#Do plant major pollinators cluster?
tip.pcols <- setNames(pcols[plants.majorpoll], 
                      sort(plants.tree$tip.label))
node.pcols<- setNames(c(tip.pcols[plants.tree$tip.label], rep("black",plants.tree$Nnode)),
                      1:(length(plants.tree$tip)+plants.tree$Nnode))

phylomorphospace(plants.tree, scent.mds.taxa, control=list(col.node=node.pcols), node.size=c(0,2), label="off", lwd=1)
text(scent.mds.taxa[,1], scent.mds.taxa[,2], rownames(scent.mds.taxa), cex=0.8, offset=0.5, xpd=T, col=tip.pcols)

library(vegan3d)
orditree3d(scent.mds.taxa, as.hclust(force.ultrametric(multi2di(plants.tree))), lwd=2, pch=19, col=tip.pcols)

phenogram(plants.tree, setNames(plants.numlinks,sort(plants.tree$tip.label)))
   Optimizing the positions of the tip labels...

png(file="ppm-%04d.png",width=600,height=600,res=120)
par(mar=c(2.1,2.1,1.1,1.1))
project.phylomorphospace(tree=plants.tree, X=scent.mds.taxa, xlab="", ylab="", node.size=c(0,0), lwd=1, direction="to", nsteps=100, fsize=0.6)
dev.off()
system("convert -delay 10 -loop 2 *.png phylochemo.gif")
file.remove(list.files(pattern=".png"))

Phylogeny scent bar chart

ccol = sample(rainbow(ncol(scent.net)))
par(mfrow=c(1,2))
par(mar=c(0.8,0,0.8,0))
plot(plants.tree, cex=1, tip.color = tip.pcols)
par(mar=c(0,0,0,0))
barplot(t(scent.net), col=ccol, names.arg=rep("",nrow(scent.net)), horiz=TRUE) 

par(mfrow=c(1,2))
par(mar=c(0.8,0,0.8,0))
plot(plants.tree, cex=1, tip.color = tip.pcols)
par(mar=c(0,0,0,0))
barplot(t(decostand(scent.net, method="tot")), col=ccol, names.arg=rep("",nrow(scent.net)), horiz=TRUE) 

Indicator compounds

library(indicspecies)
mp <- multipatt(as.data.frame(scent.net), plants.majorpoll, control=how(nperm=999))
summary(mp)
   
    Multilevel pattern analysis
    ---------------------------
   
    Association function: IndVal.g
    Significance level (alpha): 0.05
   
    Total number of species: 378
    Selected number of species: 2 
    Number of species associated to 1 group: 2 
    Number of species associated to 2 groups: 0 
    Number of species associated to 3 groups: 0 
   
    List of species associated to each combination: 
   
    Group Lepidoptera  #sps.  2 
                              stat p.value  
   un 68, 67, 81, 54, 41, 95 1.000   0.027 *
   3-carene                  0.986   0.014 *
   ---
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(paco)

# start the paco procedure
poll_phy <- polls.classtree$phylo
pla_phy <- as.phylo(hclust(vegdist(scent.net), "mcquitty"))
int <- net.lump
H <- cophenetic(polls.classtree$phylo)
P <- cophenetic(as.phylo(hclust(vegdist(scent.net), "mcquitty")))
D <- prepare_paco_data(H, P, net.lump)
   Warning in prepare_paco_data(H, P, net.lump): The HP matrix should have
   hosts in rows. It has been translated.
D <- add_pcoord(D, correction='cailliez')
# now we are ready for cophylogenetic analysis
#D <- PACo(D, nperm=1000, seed=13, method='quasiswap', symmetric=TRUE)
#save(D, file="D.Rdata")
load("D.Rdata")
# and to investigate the contribution of individual links
res <- residuals_paco(D$proc)

# to visualise this we use the ape function cophyloplot weighted by interaction contribution
# first we must make a list out of the interaction matrix
assoc <- data.frame(pol=rownames(net.lump)[which(net.lump>0, arr.ind=TRUE)[,'row']], pla=colnames(net.lump)[which(net.lump>0, arr.ind=TRUE)[,'col']])
# to weight the interactions we use the cophylogenetic contribution transformed to best show
# the differences graphically
weight <- (res^-2)/50

mnet.lump <- melt(net.lump) #replaces assoc
mnet.lump <- mnet.lump[mnet.lump$value>0,]
#cophylo with interaction strength
cophyloplot(pla_phy, poll_phy, mnet.lump[,1:2], lwd=sqrt(mnet.lump[,3])/4, show.tip.label=F,  col='steelblue', length.line=0, gap=-20, space=60)

#cophylo with residuals
cophyloplot(pla_phy, poll_phy, mnet.lump[,1:2], lwd=weight/8, show.tip.label=F, col='steelblue', length.line=0, gap=-20, space=60)

Do generalist plants have higher diversity of scent compounds?

plants.numlinks  <- rowSums(decostand(net.lump, "pa"))
polls.classtree <- compute.brlen(polls.classtree$phylo,1)
library(picante)
   Loading required package: nlme
   
   Attaching package: 'nlme'
   The following object is masked from 'package:sna':
   
       gapply
plants.pdlinks   <- pd(net.lump, polls.classtree)$PD
plants.numcompounds <- rowSums(decostand(scent.net, "pa"))
plants.sumcompounds <- rowSums(scent.net)
plants.shannoncompounds <- diversity(scent.net, index="shannon")
plants.shannonlinks <- diversity(net.lump, index="shannon")

plot(plants.numcompounds~plants.numlinks)

plot(plants.numcompounds~plants.pdlinks)

plot(plants.shannoncompounds~plants.pdlinks)

plot(plants.shannoncompounds~plants.shannonlinks)

plot(sqrt(plants.sumcompounds)~plants.pdlinks)

plot(plants.shannoncompounds~plants.majorpoll)

plot(plants.numcompounds~plants.majorpoll)

plot(sqrt(plants.sumcompounds)~plants.majorpoll)

qplot(plants.pdlinks, plants.numcompounds) + geom_smooth(method="lm")

summary(lm(plants.numcompounds[-9]~plants.pdlinks[-9]))
   
   Call:
   lm(formula = plants.numcompounds[-9] ~ plants.pdlinks[-9])
   
   Residuals:
       Min      1Q  Median      3Q     Max 
   -19.665 -10.691  -1.845   9.243  30.240 
   
   Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
   (Intercept)        20.76015    4.15181   5.000 1.33e-05 ***
   plants.pdlinks[-9]  0.02646    0.07241   0.365    0.717    
   ---
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
   
   Residual standard error: 13.82 on 38 degrees of freedom
   Multiple R-squared:  0.003502,   Adjusted R-squared:  -0.02272 
   F-statistic: 0.1335 on 1 and 38 DF,  p-value: 0.7168
data(phylocom)
pd(phylocom$sample, phylocom$phylo)
           PD SR
   clump1  16  8
   clump2a 17  8
   clump2b 18  8
   clump4  22  8
   even    30  8
   random  27  8